#include <iostream>
#include <fstream>
#include <sstream>
#include <vector>

#include <stdio.h>
#include <time.h>
#include <cmath>
#include <string>
#include "Eigen/Dense"
#include <cmath>

using namespace std;
using namespace Eigen;

#define PI 3.1415926535897932384626433832795
#define RunTime1 15e5
#define RunTime2 2e6
#define Level 57
#define NoC 56  //(level-1)
#define NoM 112 // 2*(level-1)
#define NoC1 48 //(level-1)
#define NoM1 96 // 2*(level-1)
#define freq 50.0
#define M 6
#define N 7
#define Ma 14
#define Na 15
#define Mb 2
#define Nb 3

#define NoSTN 2
#define NLine 1

using namespace std;

class IGBTsimpolorer
{
public:
    double G[8][8] = {0};
    double f[8] = {0};
    double dt = 1e-8;
    double V_nodeC, V_node1, V_nodeD, V_nodeG, V_nodeS, V_nodeAUX, V_nodemidFWD,
        V_nodeE, h = dt, Vs, VGS_old, QCE_old, QCG_old, Qeb_old,
                 Itotal, VGS, QCE, QCG, Qeb;
    void solve();
};

double max(vector<double> &data)
{
    int max_index = 0;
    for (int i = 1; i < data.size(); i++)
    {
        if (data[max_index] < data[i])
            max_index = i;
    }
    return data[max_index];
}

int main()
{
    ofstream fout("./output.csv", ios::out);
    fout << "V_nodeC,V_node1, V_nodeD, V_nodeG, V_nodeS, V_nodeAUX, V_nodemidFWD,V_nodeE,Itotal,I_my,breaktype,count" << endl;

    IGBTsimpolorer igbt;
    double R = 1e-6;
    double G = 1.0 / R;
    double r = 1002;
    double g = 1.0 / r; // 电源内阻
    double Vsource = 1.23e3;
    double Isource = Vsource *g;

    Matrix<double, 8, 8> Y;
    Matrix<double, 8, 1> J;
    Matrix<double, 8, 1> V, V_old;

    Y.setZero();
    J.setZero();
    V.setZero();
    V_old.setZero();

    vector<Matrix<double, 8, 1>> V_csv;

    for (double simulation_time = 0; simulation_time < 4000 * igbt.dt; simulation_time += igbt.dt)
    {
        igbt.VGS_old = igbt.VGS;
        igbt.QCE_old = igbt.QCE;
        igbt.QCG_old = igbt.QCG;
        igbt.Qeb_old = igbt.Qeb;

        vector<double> err;
        int count = 0;
        int breaktype = 0;
        double max_err;
        static int count1 = 0;
        static int sig = 1;
        if (count1 >= 400)
        {
            count1 = 0;
            sig = 1 - sig;
        }
        igbt.Vs = sig * 15 + (-10) * (1 - sig);
        igbt.Vs = 15;
        count1++;

        while (1)
        {
            igbt.solve();
            for (int i = 0; i < 7; i++)
            {
                J(i + 1) = igbt.f[i];
                for (int j = 0; j < 7; j++)
                {
                    Y(i + 1, j + 1) = igbt.G[i][j];
                }
            }
            Y(0, 0) += g + G;
            Y(0, 1) += -G;
            Y(1, 1) += G;
            Y(1, 0) += -G;

            J(0) += Isource;

            V = Y.inverse() * J;

            igbt.V_nodeC = V(1);
            igbt.V_node1 = V(2);
            igbt.V_nodeD = V(3);
            igbt.V_nodeG = V(4);
            igbt.V_nodeS = V(5);
            igbt.V_nodeAUX = V(6);
            igbt.V_nodemidFWD = V(7);
            igbt.V_nodeE = 0;

            if (simulation_time < igbt.dt)
            {
                breaktype = 0;
                break;
            }

            count++;
            if (count >= 1e6)
            {
                breaktype = 1;
                break;
            }

            for (int i = 0; i < 8; i++)
            {
                err.push_back(fabs(V(i) - V_old(i)));
            }
            max_err = max(err);
            // cout << "V: " << V.transpose() << endl;
            // cout << "V_old: " << V_old.transpose() << endl;
            if (max(err) < 1.0e-6)
            {
                breaktype = 2;
                break;
            }

            V_old = V;
            err.clear();

            // V_csv.push_back(V);
        }
        for (int i = 0; i < 7; i++)
        {
            fout << V(i + 1) << ",";
        }
        fout << igbt.V_nodeE << ",";
        fout << igbt.Itotal << ",";
        fout << -V(0) * g + Isource << ",";
        fout << breaktype << ",";
        fout << count << endl;
    }
    // for (int i = 0; i < V_csv.size(); i++)
    // {
    //     fout << V_csv[i].transpose() << endl;
    // }
}

void IGBTsimpolorer::solve()
{
    // double VP, K, KLM, A_FET, M_FET, N_FET, BN, M_BJT, ISAT_BJT, Vsat, Isat, VDS, ID;
    // double GmosVGS, GmosVDS, Imoseq, ID_Factor1, ID_Factor2, ID_Factor3, d1VDS, d2VDS, d3VDS, d2VGS, d3VGS;
    // double k, q, T, VT, VEB, IB, Gbjt, IBeq, RP_BJT, M_FWD, ISAT_FWD, IFWD, GFWD, IFWDeq, VFVD;
    // double C_GE, GCGS, ICGS_eq, Rg;
    // VP = 6;
    // K = 10;
    // KLM = 0;
    // A_FET = 0.4;
    // M_FET = 1.2;
    // N_FET = 1.7;
    // BN = 15;
    // M_BJT = 2;
    // ISAT_BJT = 1E-9;

    VGS = V_nodeG - V_nodeS;

    // 2.2 MOSFET
    Vsat = A_FET * pow((VGS - VP), M_FET);
    Isat = K / 2 * pow((VGS - VP), N_FET);
    VDS = V_nodeD - V_nodeS;

    if (VDS < 0)
    {
        VDS = 0;
    }

    if (VGS <= VP)
    {
        ID = 0;

        GmosVGS = 0;
        GmosVDS = 0;
        Imoseq = 0;
    }
    else if (VDS < Vsat)
    {
        ID_Factor1 = 1 + KLM * VDS;
        ID_Factor2 = 2 - VDS / Vsat;
        ID_Factor3 = VDS / Vsat;
        ID = Isat * ID_Factor1 * ID_Factor2 * ID_Factor3;

        d1VDS = KLM;
        d2VDS = -1 / Vsat;
        d3VDS = 1 / Vsat;
        GmosVDS = Isat * (d1VDS * ID_Factor2 * ID_Factor3 + ID_Factor1 * d2VDS * ID_Factor3 + ID_Factor1 * ID_Factor2 * d3VDS);
        d2VGS = VDS / pow(Vsat, 2) * A_FET * M_FET * pow((VGS - VP), (M_FET - 1));
        d3VGS = -VDS / pow(Vsat, 2) * A_FET * M_FET * pow((VGS - VP), (M_FET - 1));
        GmosVGS = K / 2 * N_FET * pow((VGS - VP), (N_FET - 1)) * ID_Factor1 * ID_Factor2 * ID_Factor3 + Isat * ID_Factor1 * d2VGS * ID_Factor3 + Isat * ID_Factor1 * ID_Factor2 * d3VGS;
        Imoseq = ID - GmosVDS * VDS - GmosVGS * (VGS);
    }
    else
    {
        ID = Isat * (1 + KLM * VDS);

        GmosVDS = Isat * KLM;
        GmosVGS = K / 2 * N_FET * pow((VGS - VP), (N_FET - 1)) * (1 + KLM * VDS);
        Imoseq = ID - GmosVDS * VDS - GmosVGS * (VGS);
    }

    // 2.3 BJT
    // k = 1.381E-23;
    // q = 1.602E-19;
    // T = 300;
    // VT = k * T / q;
    VEB = V_node1 - V_nodeD;

    if (VEB > 2)
    {
        VEB = 2;
    }

    IB = ISAT_BJT * (exp(VEB / M_BJT / VT) - 1);
    Gbjt = ISAT_BJT * exp(VEB / M_BJT / VT) / M_BJT / VT;
    IBeq = IB - Gbjt * VEB;

    RP_BJT = 125;

    ////////////2.4 FWD

    // M_FWD = 1.21855;
    // ISAT_FWD = 1.001e-6;
    VFVD = V_nodeAUX - V_nodemidFWD;
    ///////////////////////////////////important////////////////////////////////////
    if (VFVD > 1)
    {
        VFVD = 1;
    }
    ////////////////////////////////////////////////////////////////////////////////
    IFWD = ISAT_FWD * (exp(VFVD / M_FWD / VT) - 1);
    GFWD = ISAT_FWD * exp(VFVD / M_FWD / VT) / M_FWD / VT;
    IFWDeq = IFWD - GFWD * VFVD;

    // double RB_FWD = 0.0001;
    // double INOM = 3000;
    double R_FWD = RB_FWD / sqrt(1 + IFWD / INOM);
    double G_FWD = 1 / R_FWD;

    /////////////2.5 C_GE
    // C_GE = 0.7e-7;
    // GCGS = C_GE / h;
    ICGS_eq = -GCGS * VGS_old;

    /////////////2.6 C_CE
    // double C0_CE = 2e-8;
    // double VDIFF_CE = 0.2;
    // double VSHIFT_CE = 0;
    // double ALPHA_CE = 0.5;
    // double BETA_CE = 2;
    // double DELTA_CE = 0.0001;

    // double VCE, B, VVJNCT_CE, CEa, CEb, AA_CE, BB_CE;

    VCE = V_node1 - V_nodeS; // QCE_old
    B = VSHIFT_CE * VDIFF_CE;
    VVJNCT_CE = B - VCE;
    CEa = BETA_CE - 1;
    CEb = 1 - DELTA_CE;

    AA_CE = (ALPHA_CE * CEb / CEa) / VDIFF_CE;
    BB_CE = (1 - DELTA_CE) * VDIFF_CE / (1 - ALPHA_CE);

    // double CCE, ICE, GCE, ICEeq;
    if (VVJNCT_CE >= 0)
    {
        CCE = C0_CE * (1 + CEa * (1 - exp(-VVJNCT_CE * AA_CE)));

        // Charge Calculation by Integration of CCE(from 0 to VCE)
        QCE = C0_CE * (VCE + CEa * (VCE - exp(VCE * AA_CE) / AA_CE)) + C0_CE * CEa / AA_CE;
        ICE = (QCE - QCE_old) / h;
        GCE = CCE / h;
        ICEeq = ICE - GCE * VCE;
    }

    else
    {
        CCE = C0_CE * (DELTA_CE + CEb / pow((1 - VVJNCT_CE / VDIFF_CE), ALPHA_CE));
        // Charge Calculation by Integration of CCE(from 0 to VCE)
        QCE = C0_CE * (DELTA_CE * VCE + BB_CE * pow((1 + VCE / VDIFF_CE), (1 - ALPHA_CE))) - C0_CE * BB_CE;
        ICE = (QCE - QCE_old) / h;
        GCE = CCE / h;
        ICEeq = ICE - GCE * VCE;
    }
    ///////////////////2.7 C_CG
    // double C0_CG = 5e-8;
    // double VDIFF_CG_OFF = 0.1;
    // double VDIFF_CG_ON = 0.1;
    // double VSHIFT_CG_OFF = 0;
    // double VSHIFT_CG_ON = 0;
    // double ALPHA_CG_OFF = 0.5;
    // double ALPHA_CG_ON = 0.5;
    // double BETA_CG_OFF = 2;
    // double BETA_CG_ON = 2;
    // double DELTA_CG_OFF = 0.0001;
    // double DELTA_CG_ON = 0.0001;

    // double VCG, VVJNCT_CG, CGa, CGb, AA_CG, BB_CG, CCG, ICG, GCG, ICGeq;

    VCG = V_node1 - V_nodeG; // QCG_old
    if (VGS <= VP)
    {
        B = VSHIFT_CG_OFF * VDIFF_CG_OFF;
        VVJNCT_CG = B - VCG;
        CGa = BETA_CG_OFF - 1;
        CGb = 1 - DELTA_CG_OFF;

        AA_CG = ALPHA_CG_OFF * CGb / CGa / VDIFF_CG_OFF;
        BB_CG = (1 - DELTA_CG_OFF) * VDIFF_CG_OFF / (1 - ALPHA_CG_OFF);

        if (VVJNCT_CG >= 0)
        {
            CCG = C0_CG * (1 + CGa * (1 - exp(-VVJNCT_CG * AA_CG)));
            // Charge Calculation by Integration of CCE(from 0 to VCE)
            QCG = C0_CG * (VCG + CGa * (VCG - exp(VCG * AA_CG) / AA_CG)) + C0_CG * CGa / AA_CG;
            ICG = (QCG - QCG_old) / h;
            GCG = CCG / h;
            ICGeq = ICG - GCG * VCG;
        }
        else
        {
            CCG = C0_CG * (DELTA_CG_OFF + CGb / pow((1 - VVJNCT_CG / VDIFF_CG_OFF), ALPHA_CG_OFF));
            // Charge Calculation by Integration of CCE(from 0 to VCE)
            QCG = C0_CG * (DELTA_CG_OFF * VCG + BB_CG * pow((1 + VCG / VDIFF_CG_OFF), (1 - ALPHA_CG_OFF))) - C0_CG * BB_CG;
            ICG = (QCG - QCG_old) / h;
            GCG = CCG / h;
            ICGeq = ICG - GCG * VCG;
        }
    }
    else
    {
        B = VSHIFT_CG_ON * VDIFF_CG_ON;
        VVJNCT_CG = B - VCG;
        CGa = BETA_CG_ON - 1;
        CGb = 1 - DELTA_CG_ON;

        AA_CG = ALPHA_CG_ON * CGb / CGa / VDIFF_CG_ON;
        BB_CG = (1 - DELTA_CG_ON) * VDIFF_CG_ON / (1 - ALPHA_CG_ON);

        if (VVJNCT_CG >= 0)
        {
            CCG = C0_CG * (1 + CGa * (1 - exp(-VVJNCT_CG * AA_CG)));
            // Charge Calculation by Integration of CCE(from 0 to VCE)
            QCG = C0_CG * (VCG + CGa * (VCG - exp(VCG * AA_CG) / AA_CG)) + C0_CG * CGa / AA_CG;
            ICG = (QCG - QCG_old) / h;
            GCG = CCG / h;
            ICGeq = ICG - GCG * VCG;
        }
        else
        {
            CCG = C0_CG * (DELTA_CG_ON + CGb / pow((1 - VVJNCT_CG / VDIFF_CG_ON), ALPHA_CG_ON));
            // Charge Calculation by Integration of CCE(from 0 to VCE)
            QCG = C0_CG * (DELTA_CG_ON * VCG + BB_CG * pow((1 + VCG / VDIFF_CG_ON), (1 - ALPHA_CG_ON))) - C0_CG * BB_CG;
            ICG = (QCG - QCG_old) / h;
            GCG = CCG / h;
            ICGeq = ICG - GCG * VCG;
        }
    }

    //////////////2.8 Ceb(diff)
    // double TAU_BE = 5e-8;
    // double Ceb, Ieb, GCeb, ICebeq;

    Ceb = TAU_BE * (IB + ISAT_BJT) / M_BJT / VT;
    Qeb = TAU_BE * ISAT_BJT * (exp(VEB / M_BJT / VT) - 1);
    Ieb = (Qeb - Qeb_old) / h;
    GCeb = Ceb / h;
    ICebeq = Ieb - GCeb * VEB;

    /////////////2.6 Connector

    Rg = 10.2;
    double LG = 0;
    double LC = 5e-9;
    double RC = 1.15625e-5;
    double LE = 5e-9;
    double RE = 1.15625e-5;
    double RAUX = 1.15625e-5;
    double LAUX = 1.75e-10;

    f[0] = 0;                                         // node C
    f[1] = -(BN + 1) * IBeq - ICEeq - ICGeq - ICebeq; // node 1
    f[2] = -Imoseq + IBeq + ICebeq;                   // node D
    f[3] = Vs / Rg - ICGS_eq + ICGeq;                 // node G
    f[4] = Imoseq + BN * IBeq + ICGS_eq + ICEeq;      // node S
    f[5] = -IFWDeq - Vs / Rg;                         // node AUX
    f[6] = IFWDeq;                                    // node mid FWD
    f[7] = 0;                                         // node E

    G[0][0] = 1 / RC;
    G[0][1] = -1 / RC;
    G[0][2] = 0;
    G[0][3] = 0;
    G[0][4] = 0;
    G[0][5] = 0;
    G[0][6] = 0;
    G[0][7] = 0;

    G[1][0] = -1 / RC;
    G[1][1] = (BN + 1) * Gbjt + 1 / RP_BJT + G_FWD + 1 / RC + GCE + GCG + GCeb;
    // cout << "===============" << endl;
    // printf("%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf", BN, Gbjt, RP_BJT, G_FWD, RC, GCE, GCG, GCeb);
    G[1][2] = -(BN + 1) * Gbjt - 1 / RP_BJT - GCeb;
    G[1][3] = -GCG;
    G[1][4] = -GCE;
    G[1][5] = 0;
    G[1][6] = -G_FWD;
    G[1][7] = 0;

    G[2][0] = 0;
    G[2][1] = -Gbjt - 1 / RP_BJT - GCeb;
    G[2][2] = GmosVDS + Gbjt + 1 / RP_BJT + GCeb;
    G[2][3] = GmosVGS;
    G[2][4] = -GmosVDS - GmosVGS;
    G[2][5] = 0;
    G[2][6] = 0;
    G[2][7] = 0;

    G[3][0] = 0;
    G[3][1] = -GCG;
    G[3][2] = 0;
    G[3][3] = 1 / Rg + GCGS + GCG;
    G[3][4] = -1 / Rg - GCGS;
    G[3][5] = 0;
    G[3][6] = 0;
    G[3][7] = 0;

    G[4][0] = 0;
    G[4][1] = -BN * Gbjt - GCE;
    G[4][2] = -GmosVDS + BN * Gbjt;
    G[4][3] = -GmosVGS - GCGS;
    G[4][4] = GmosVDS + GmosVGS + GCGS + 1 / RAUX + GCE;
    G[4][5] = -1 / RAUX;
    G[4][6] = 0;
    G[4][7] = 0;

    G[5][0] = 0;
    G[5][1] = 0;
    G[5][2] = 0;
    G[5][3] = -1 / Rg;
    G[5][4] = -1 / RAUX;
    G[5][5] = GFWD + 1 / Rg + 1 / RAUX + 1 / RE;
    G[5][6] = -GFWD;
    G[5][7] = -1 / RE;

    G[6][0] = 0;
    G[6][1] = -G_FWD;
    G[6][2] = 0;
    G[6][3] = 0;
    G[6][4] = 0;
    G[6][5] = -GFWD;
    G[6][6] = G_FWD + GFWD;
    G[6][7] = 0;

    G[7][0] = 0;
    G[7][1] = 0;
    G[7][2] = 0;
    G[7][3] = 0;
    G[7][4] = 0;
    G[7][5] = -1 / RE;
    G[7][6] = 0;
    G[7][7] = 1 / RE;

    Itotal = (V_nodeC - V_node1) / RC;
}